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ABSTRACT 



<^ Markov processes are shown to be consistent with metastable states seen in pulsar phe- 

nomena, including intensity nulling, pulse-shape mode changes, subpulse drift rates, spindown 

CIh rates, and X-ray emission, based on the typically broad and monotonic distributions of state 

"^ lifetimes. Markovianity implies a nonlinear magnetospheric system in which state changes 

occur stochastically, corresponding to transitions between local minima in an effective po- 
tential. State durations (though not transition times) are thus largely decoupled from the 
characteristic time scales of various magnetospheric processes. Dyadic states are common 
but some objects show at least four states with some transitions forbidden. Another case is 
the long-term intermittent pulsar B1931-I-24 that has binary radio-emission and torque states 
with wide, but non-monotonic duration distributions. It also shows a quasi-period of 38 it 5 
days in a 13-yr time sequence, suggesting stochastic resonance in a Markov system with a 



o 

?H forcing function that could be strictly periodic or quasi-periodic. Nonlinear phenomena are 



^ associated with time-dependent activity in the acceleration region near each magnetic polar 

'— ' cap. The polar-cap diode is altered by feedback from the outer magnetosphere and by re- 

j—{ turn currents from an equatorial disk that may also cause the neutron star to episodically 

^ charge and discharge. Orbital perturbations in the disk provide a natural periodicity for the 

^^ forcing function in the stochastic resonance interpretation of B1931-I-24. Disk dynamics may 

Q^ introduce additional time scales in observed phenomena. Future work can test the Markov 

interpretation, identify which pulsar types have a propensity for state changes, and clarify the 

role of selection effects. 



Introduction 



X 

^ Ever since their discovery it has been well known that radio pulses from pulsars vary dramatically on a wide 
range of time scales. By comparison, incoherent optical to gamma-ray emission shows little variability, 
apart from the basic pulsation. Consequently, radio variability was widely considered to result from 
changes in radio coherence that did not couple to high-energy emission or to the overall loss of rotational 
energy. That has now changed with the identification of large changes in spindown rates on time scales 
~ 10^ s (Kramer et al. 2006; Camilo et al.|2012 Lorimer et al.|2012) that correlate with distinct on-and-off 



radio emission states. Equally important is the detection of X-ray emission state changes on time scales 
'^ 10^-10^ s that are linked to specific radio modes in the pulsar B0943-I-10 (Hermsen et al.|2013). These 



phenomena indicate that radio emission, though energetically small, traces the dominant energy channels 
of the magnetosphere, including acceleration regions and the large scale structure of the magnetosphere 
itself. 



The ansatz for this paper is consistent with similar ideas expressed by other authors (Ruderman &; 



Sutherland 1975 Cordes 1979 Jones 1982 Cordes 1983 Jones 1986 Timokhin 2010 Jones 2011 Li et 



al. 2012c van Leeuwen & Timokhin 2012), that the observed state changes on both short and long time 



scales reflect changes in (1) the global state and energetics of a pulsar's magnetosphere, as implied by the 
large changes in spindown rate seen in intermittent pulsars; (2) the accelerating potential directly above 
the magnetic polar caps that generates relativistic particles, as evidenced by mode changes; and (3) the 
temperature of the hot polar cap due to particle bombardment as demonstrated by X-ray low and high 
states that correlate with mode changes. 

This paper addresses issues in understanding pulsar magnetospheres using the statistical framework of 
Markov chains: (1) What is an accurate yet parsimonious description of state changes in pulse sequences? 
(2) Are some or even all observed state changes connected through some common paradigm? In particular, 
do intensity nulls occurring on time scales of hours or less have the same physical origin as the much 
longer intermittency on days to week-like time scales? (3) How directly related are observed phenomena 
to physical changes in the pulsar magnetosphere? Are there underlying state changes that are not directly 
observable? (4) What determines the occurrence of a state change? Is it triggered externally or is it a 
threshold effect in the internal dynamics of the magnetosphere? 

In Section [2] we discuss observations of state changes that are relevant to our modeling. In Sections [3]J5] 
we summarize general properties of Markov processes and apply them to some of the key observations. 
In Sections |6] and [7] we discuss spin variations expected from a two-state Markov process with reference 
to the medium-term intermittent pulsar B0823-I-26 and the long-term intermittent pulsar B1931-I-24. In 
Section [8] we relate our findings to pulsar magnetospheres and their surroundings and in Section [9] we 
summarize our results and our conclusions. 



2. State Changes in Pulsars 

Table [l] itemizes some of the phenomena that we discuss, all of which have been seen only in relatively old, 
canonical pulsars but not in young pulsars, in millisecond pulsars (MSPs), or magnetars. This suggests 
that state changes occur in objects that are in a particular region of magnetic field - spin-period space, 
though it is possible that selection effects, such as the inability to study single pulses from most MSPs, 
prohibits recognition of state changes. This situation will change as more sensitive measurements are 
obtained. 

The distribution of state durations is one of the primary tools used in this paper for characterizing state 
changes. Figure [T] presents examples for six pulsars that show nulls and bursts or pulse-shape mode 
changes. Within counting errors, the histograms decrease monotonically as expected for a two-state 
Markov process (viz. Equationfl], as discussed further below). In several cases (B0834-I-06, the abnormal 
mode for B0329-I-54, and bursts for Bl 133-1-16) there is an excess of short-duration states, as also pointed 
out for other objects ( van Leeuwen et al.| 2002 Kloumann & Rankin 2010[ [Gajjar et al. 2012). The 
duration histogram for bursts is shown for one of those, B1944-I-17, in Figure^ (left panel). Apparent 
counterexamples to these trends are shown in the right-hand panel of Figure [2j which gives duration 



histograms of bursts (B) and quiet (Q) modes from B2303-I-30. However, as pointed out by Redman 



et al. (2005), whose Figure 4 is the basis for Figure 2 short-duration B and Q sequences are probably 



undercounted, causing some of the departure from an exponential-like distribution. 

While an excess of short-duration states appears to be intrinsic to pulsar emission in some cases, in others 
there is inadequate signal-to-noise ratio (S/N) to avoid false transitions between nulls and bursts caused 



by additive noise. Since noise is independent between pulse periods, the durations of such mis-identified 
states are typically only one period long. This is discussed further in Section [3j 
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Fig. 1. — Histograms of durations of nulls, bursts, and modes for six pulsars. Frames (a) and (b) show 



nuU and burst histograms for B1133+16 at 325 MHz from |Bhat etaL] ( p007| ); (c) Nulls for B0809+74 at 
350 MHz dvan Leeuwen et al.|[2002l ); (d) Nulls for B0818-13 ( jJanssen fc van Leeuweill[2004| ; (e) Bursts 
for B0834+06 ( [Redman fc Rankin] [20091 ) ; (0 Nulls for B2303+30 ( [Redman et "aLJ [2005] ) ; (g) and (h) 
Normal and abnormal mode durations for B0329+54 ( Chen et al.|2011 ). The smooth lines are Equation [l] 
evaluated using the mean state duration and the number of counts in each case; they are not least-squares 
fits and are given solely as fiducial examples. The maximum number of counts is given in each frame. 



Markov and Stochastic Resonance Models 



The broad, exponential-like distributions for state changes, along with the lack of any correlation between 
adjacent null and burst durations, are naturally explained by pulse sequences that conform to an under- 
lying Markov chain. So too is the inconsistency of nulls and bursts with the Wald-Wolfowitz 'runs' test 
(Redman & Rankin 2009). An additional phenomenon — stochastic resonance — is suggested by the 



binary state switching in the long-term intermittent pulsar B1931+24. The durations of on-and-off states 
for this object also have wide distributions but they are neither exponential nor monotonic. In addition, 
an underlying period of 38 it 5 days is discernible in a 13-yr data set ( Young et al.[[2013 ). While sampling 
incompleteness may distort the distributions, the phenomena are indicative of stochastic resonance in a 
nonlinear system that has binary states and is driven by a forcing function that could be strictly periodic 



Table 1. State Change Phenomena in Pulsars 



Phenomenon 



Properties 



Number of 
States 



Examples 



References 



Nulls and bursts 



On:off flux density ratio > 100. 
Null fractions: /„ < 94%. 
Durations < 1 s to hours. 
Residence time PDFs: 

1. broad, monotonic 

2. dual components (often). 



2 to 4 



B0031-07 
B1237+25 
B1944+17 
B0826-34 



4,5 



Mode changes 



Changes in average pulse shapes. 



2 to 3 



B0329+54 
B1237+25 



6,7 



Subpulse drifts 



Systematic pulse-phase drifts 

through pulse window. 
Preferred drift rates. 
Allowed, forbidden state sequences. 
Mixed with unorganized drift states. 
Quasi-periodicities in > 55% of pulsars. 



<4 



B0031-07, 
B1822-09 
B1944+17 
B2319+60 



10, 11 



Long-term intermittency 

Profile and torque changes 
Radio/X-ray correlation 
Extreme bursts 



Off states: hours to years. 2 

Torque 1.5 to 2.5 times larger in on state. 

38-day quasi-periodicity in B1931+24 

Correlated mode and torque changes. > 2 

Month-like time scales. 

X-ray on-off state transitions 2 

correlate with radio modes. 

Giant pulses. 2 

Sporadic bursts with underlying period. 2 



J1832+0029 
J1841-0500 
B1931-h24 



B0943-(-10 



12 
13 
14, 15 



Crab pulsar 18 

(J0534+2200) 

Rotating radio 19,20,21,22 

transients 



References. — 1. |Ritchings|([T976|; 2. |Biggs| ([T992| ; 3. |Deich et al.|([T986|; 4. [Kloumann fc Rankm|(|2010|; 5. [Gajjar 



et al. 



20121; 6 



10. Backer 



al. 



Bartel et al. 



(19821; 7. Wang et al. 



(19731; 11. Weltevrede et al. (20061; 12 



Lorimer et al. 



(20121; 13 



McLaughlin et al.| ( |2006| ); 20. [Deneva et al.| ( |2009[ ); 21. [Keane et al.| ( |2011[ ); 22. [Palliyaguru et al.| ( |2011 l 



20071; 8. Huguenin et al. (19701; 9. Wright fc Fowler (1981bl; 



Camilo et al. 



(20121; 14 



Kramer et 



( |2006| ); 15. [Young ct al.| ( |2013[ ); 16. |Lyne etaT] ( |2010| ); 17. [Hermsen et al] ( |2013| ); 18. [Lundgren et al.| ( |1995| ); 19 
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Fig. 2. — Left: Histogram of burst durations for B1944+17 at 327 MHz (after Kloumann & Rankin 
2010). The solid line shows the predicted distribution for a four-state Markov model where two burst 
states are counted together as discussed in the text (Section |4.2[). Right: Histograms of the durations 



of the burst (B) and quiet (Q) modes for B2303-I-30 from Redman et al. (2005). The smooth lines are 
Equation [T] evaluated using the mean state duration and the number of counts in each case. 



or itself quasi-periodic. 

Along with binary states, some pulsars show multiple discrete values of subpulse drift rate or quasi-periods 
in pulse amplitudes. These have been compellingly associated with a carousel of sub-beams that rotate 



through the overall radio emission cone with circulation times of tens to hundreds of seconds ( Deshpande 



&: Rankiii||1999 Rankin &: Stappers||2008 ). The changes in drift rate or quasi-period, often combined with 
nulling, imply that some objects display at least four states. These are also easily describable as Markov 
chains but with some state transitions disallowed, as discussed below. 

Our approach is to generally discuss observed state-switching in terms of Markov processes and stochastic 
resonance. 



3.1. Markov Processes 



Markov processes are summarized here with reference to Papoulis (1991). We consider n-state Markov 
processes (denoted M„) that have state values S = {si,S2,--- ,5^}. Transitions between states are 
described by an n x n stochastic matrix Q whose elements are the probabilities qij,i,j = 1, . . . ,n for 
changing from state Sj to state Sj in a single time step. For observations we discuss, the time step is a 
single rotation periodr Normalization across a row is ^ ■ qij = 1 and the probability for staying in the 



■th 



state is the "metastability" qu. The duration T of the i state is a random variable with probability 
(X qj^ for T time steps, mean Tj = (1 — qu)'^, and rms ctt^ = Ti^/qu. The probability density function 



(PDF) of T for the i state is given by g/j = (1 



TfY 



normalized by the sum of this quantity over all 



^Ultimately one must consider variations that occur on shorter time scales in the corotating frame and include the 
stroboscopic effect of pulse-window sampling related to rotation. This kind of analysis is deferred to another paper. 



values of r, or 



fT,{T)=Tr^{l-TrY \ T = l,2,---. (1) 



The PDF has constant slope dln/y. (r)/(ir = ln(l — 1/Tj) and is generally steeper than an exponential 
that has the same mean, Tj. However, as Tj — )• cxd (or equivalently qu — )• 1), the PDF tends to a one-sided 
exponential function with 0"^. = Tj. 

The ensemble probabilities of each state are given by the state probability vector P = {pi,P2,- ' ' ,Pn)- 
Stationarity, which we assume, requires that PQ = P, so P is the left eigenvector of Q with unit 
eigenvalue. The transition matrix after t time steps, Q*, converges to a form Q°° whose rows are all equal 
to the state probability vector. The convergence time is some multiple of the longest-duration state. 

A two-state M2 process with S = {n, b} for nulls and bursts, for example, has pi^2 = ^1,2/(^1 + T2). The 
nulling fraction equals pi and can be written as 

f — ^1 _ 1 ~ 922 _ g21 ^„^ 

T1+T2 2 - qn - q22 qi2 + 921 

making explicit that /„ can be identical for objects that have markedly different mean state durations, 
as observed. Examples of two-state Markov processes in Figure p^ (top three rows) show the dependence 
of the time series and state-duration histograms on the metastabilities qu, ^22- 

To model a time series, the number of states and the transition matrix need to be determined. There 
are n{n — 1) unique elements in Q for an n-state process, after accounting for normalization across rows. 
Of these, n elements can be determined from the mean state durations, Tj,i = 1, . . . ,n. Frequencies of 
occurrence of the states constrain an additional n — 1 elements. The remaining n{n — 3) + 1 elements need 
to be obtained using additional assumptions or measurements. For n = 2, the system can be determined 
solely from the mean state durations so measured frequencies of occurrence provide redundancy. For 
n = 3, mean durations and frequencies need to be combined with one additional measurement is needed. 
We will also consider n = 4 processes that require five additional elements. 



3.2. Nonlinear Systems and Stochastic Resonance 

Markov processes are often used to describe non-linear systems that have an underlying potential-energy 
function with minima corresponding to the Markov states (e.g. Becker &: Rein ten Wolde||2012). Random 



motions ('noise') within potential wells induce state changes at rates that depend on well depths and 
noise strengths. State changes can also be driven by a forcing function that itself is stochastic, chaotic, 
or deterministic (Anishchenko, Anufrieva &: Vadivasova||2006 ) . A strictly periodic forcing function acting 



in concert with noise can induce state changes that are only quasi-periodic (Gammaitoni et al. 1998), a 
phenomenon known as stochastic resonance. 

When driven by both noise and a forcing function d{t), a system potential having two local minima 
separated by a barrier can be described as an M2 process in the adiabatic regime where d(t) is slowly 
varying. The transition matrix is 

^ _ / 1 - <Zi2e^^'^W q,,e^^m \ 

where Ai^2 ^ depend on properties of the wells and noise. Generally, the process is non-stationary. 
For A12 = 0, Q reverts to the form for a simple M2 process with time independent elements. Later we 



consider periodic forcing with d{t) = cos(27rt/Pf + </>) where Pf is the period and (j) is an arbitrary phase. 
Over a cycle of d{t), the well depths oscillate with the two wells out of phase by vr. 

Examples of two-state Markov processes in Figure[3]with different degrees of stochastic resonance (bottom 
three rows) show how state-duration histograms depend on the system parameters. Unlike the top three 
rows that show exponential-like histograms, stochastic resonance makes the histograms non-monotonic 
and periodic. Power spectra shown in the right-hand panels are featureless in the top three rows but show 
a spectral line to varying degrees in the bottom three rows at a frequency ~ 20 mcycles step~^. 
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Fig. 3. — Time series, histograms of state durations, and power spectra for two-state processes. Left: 
single-realization time series for different values of the metastability probabilities qu and g22 and for the 
stochastic resonance amplitudes vli,2) as given in the top of each frame. The top three rows are for pure 
Markov processes {Ai^2 = 0) and the bottom three include a periodic forcing function with Ai 2 > and 
period Pf = 50 time steps. Center: Histograms of the durations of low states and high states states based 
on 5000 realizations of the time series. The logarithmic vertical scales are from 1 to 10^'^. Right: Power 
spectra averaged over 5000 realizations plotted up to one-half the Nyquist frequency. The logarithmic 
vertical scale ranges from 10~^ to 0.1. The 50-step period for the stochastic resonance in the lower three 
panels is evident to varying degrees as a spectral line at a frequency of 20 milli-cycles step~^. 



4. Application to Short-term Nulling and Pulse-Profile Mode Changes 

Nulls and bursts are often described empirically as two-state phenomena and so are some pulse-shape 
'mode' changes. The null fraction appears to be characteristic of a pulsar and is constant over decades in 
well-studied objects, suggesting that the transition matrix in a Markov model is stable in the same manner 
as the average pulse profile. Histograms of state durations are almost always broad and monotonically 



decreasing (Deich et al. 1986 Redman et al. 2005 Bhat et al. 2007 Kloumann & Rankin 2010 IChen et 



ai:] [20TH [Gajjar et al.|[2()12D with many (e.g. B0835-41, B1112-h50, B1133+16, B2111+46, B2303+30) 
showing consistency with the exponential- like form expected for an M2 Markov process. 

Null and burst sequences fail the Wald-Wolfowitz 'runs' test that is predicated on transitions occurring 



independently of the current state, as in coin tossing (Redman &: Rankin 2009) or any independent, 
identically distributed (i.i.d.) process. Coin tossing corresponds to a Markov process whose transition 
matrix Q has identical rows. However, transition matrices consistent with observed nulls and bursts 
necessarily have different rows and will fail the runs test. The test statistic used for the runs test is 
Z = [R — {R))/aR (Redman &: Rankin 2009), where R is the length of a null or burst; Z has zero mean 



and unit variance for an i.i.d. process. In general an M2 process has mean (Z) = ^/N{T-^ + ^2 ~ ^) 



N [(1 — fn)Ti — 1] where A^ ^ 1 is the number of pulses analyzed. An i.i.d. process has T^ +^2^ = 1 
so that {Z) = 0. However a Markov process with nulling fraction /„ = 0.32 and A^ = 1024 pulses (an 
example in Redman & Rankin (2009)) yields {Z) = —27.3. It is noteworthy that estimates of Z for 25 



out of 26 pulsars range between -0.5 and -43 (Redman & Rankin 2009 Gajjar et al. 2012), as expected 



for an M2 process that accounts for the observed state durations. 



4.1. False Transitions from Measurement Noise and Beaming 

In this section we consider the effects of misidentified state transitions that inevitably result from mea- 
surements with finite S/N. Coherent emission from pulsars typically shows a broad pulse-amplitude dis- 
tribution (e.g. Gajjar et al.| 2012 Figure 2), often taking a log-normal form (e.g. Burke-Spolaor et al. 



2012) for the phase-integrated intensity and in some cases a power-law form (e.g. Lundgren et al. 1995 
Cognard et al]|1996t [Cordes et al.||2004D . 

A model for the averaged intensity (pulse 'energy') in an on-pulse window is 



I{t) = b{t)A{t\M2{t)) + N{t), 



(4) 



where N{t) is white, Gaussian noise with variance o"^ and A{t\M2{t)) is a state-dependent pulse am- 
plitude. The time-dependent beaming function h{t) is a multiplicative factor that models pulse-to-pulse 
modulations caused by the entire radiation beam being only partially and stochastically filled with sub- 



beams or containing periodically spaced beams, as in the carousel model of Deshpande &: Rankin (1999). 
Letting the two states be S = {0, 1}, the pulse amplitude is 



A{t\M2{t) 



M2 = 

a(/i,a)/(a) M2 = 1, 



(5) 



where a{fj,,a) is a log-normal random variable with parameters /x and a. When M2 = 1, a statistically 
independent value is drawn for a(^, a) and normalized by the mean, (a). Defined this way, {A{t)) = 1 for 
the M2 = 1 state and the average signal to noise ratio is S/N = (A)/(Tn = I/cn- 



4-1.1. False Transitions from Additive Noise 

First we consider the effects of additive noise when the beaming function is constant. Figures [4] and [5] 
show time series and histograms for simulations using Equation |4j In Figure |4] the low and high states 
are well separated so that additive noise does not affect the identification of the two states in a time 
series. Time series were generated according to Equations E]J5] using S/N = 20, fi = 1.5, and a = 0.2 and 
a threshold of x = 3o"7v was used to discriminate between low and high states. The lower S/N case in 
Figurejsjwith S/N = 10, // = 0.8, and a = 0.5 shows overlap between the two observable state amplitudes, 
with consequent misidentified transitions. These alter the duration histograms by adding many short nulls 
and shortening the durations of bursts, which respectively cause the narrow feature near T = in the 
inferred-null histogram and the steepening of the inferred burst histogram in the lower-right of the figure. 

The false-transition probabilities for a threshold x and for PDFs /Ar(A^) and /a (a) for the noise and pulse 
amplitudes, respectively, are 

oo 



r-12 = / dNfN{N) (6) 

poo px—a 

r2i = / dafaia) / dN fN{N). 



Given the statistical independence of Markov transitions and noise-induced threshold crossings, the com- 
bined off-diagonal transition probabilities are (using 'M' to denote the noise- free Markov values) 

qi2 = qi2,M + ri2 - qi2,M'^l2 > gi2,M, (7) 

q2i = q2i,M + r2i - q2i,Mr2i > q2i,M, (8) 

showing that the interstate transition probabilities are always increased by additive noise, as expected. 

For the case presented in Figure [5} r2i > ri2 because the Gaussian noise PDF falls off more rapidly for 
high values than does the log-normal PDF going to small values. This causes the spike at short durations 
to appear only in the null-duration PDF. For other choices of PDFs, S/N, and threshold, spikes can 
appear in one or the other (or both) of the null and burst PDFs . 



4-1-2. Pseudo-Nulls from the Multiplicative Beaming Function 



Kloumann Sz Rankin (2010 hereafter KLIO) identified an excess number of short nulls and short bursts 
over what would be expected from an M2 process having long mean state durations. Their physical 
explanation is that the excess short nulls are "pseudo" nulls resulting from the granularity of conal 
carousel beams while true nulls (both short and long) are caused by bona fide state changes, either a 
temporal cessation or a strong diminution of the radio intensity. This interpretation implies that the two 
processes (beaming and actual nulling) are independent and therefore multiplicative. 

When the beaming function b{t) is included in the intensity model of Equation E^ as a highly modulated 
quantity with minima well below the mean, the statistics of true nulls and bursts are modified. Short, 
pseudo nulls ~one period in duration truncate bursts and therefore decrease the mean durations of both 
nulls and bursts while increasing their numbers in a time series. It can be shown that any modulation, 
periodic or otherwise, will have the same effect on the null and burst duration histograms. It is not so 
clear that beam modulations will have the analogous effect of producing short bursts unless there are 
occasional times when b{t) is especially large for durations of one spin period (or less). 
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Fig. 4. — Simulated time series and histograms for a case where all pulse amplitudes are well above the 
threshold for dividing null and burst pulses. The left column shows a single time-series realization while 
the histograms in the right-hand column are based on 1000 realizations. Top left: time series of the 
Markov process. Middle left: time series of the intensity that includes additive noise and pulse-amplitude 
variations (Equation H|; the mean signal-to- noise ratio and mean state durations are indicated. Lower 
left: inferred state sequence obtained by identifying transitions in the intensity time series. Top right: 
histogram of pulse amplitudes; the dashed line represents the 3a threshold for separating null and burst 
pulses. Middle right: histograms of the true durations of null and bursts states; the estimated mean 
durations are indicated. Bottom right: histograms of the null and burst durations inferred from the 
intensity. 



4.2. Dual-Component Nulls and Bursts 

Several objects in addition to B1944+17 show dual characteristic durations for nulls or bursts. Some but 
not all short nulls and bursts with ~ 1 period durations could be caused by false transitions from additive 
noise or by pseudo- nulls due to multiplicative beaming variations, as discussed above. 

Here we acknowledge such contaminating effects on null and burst statistics but proceed with a Markov 
description, which is certainly allowed no matter the cause of the various kinds of nulls and bursts. In 
particular, we demonstrate that the estimated distributions for null and burst durations can be accounted 
for by an M4 process. The state set S = {ni,n2, 61, 62} comprises two null states (ni_2) and two burst 
states (61,2) that are not observationally distinguished, apart from their durationsrl To proceed, we adopt 



Models that include multiple states that are indistinguishable observationally are hidden Markov models (HMM; Rabiner 
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Fig. 5. — Same as Figure |4] but for a case where the pulse amplitude distribution and additive noise 
distributions overlap significantly. 



a minimalist model where transitions are allowed only between a null state and a burst state and vice versa 
but ni - 722 and 61-62 transitions in either direction are forbidden. Of the 16 equations needed to solve for 
the matrix elements, four come from normalization across rows, four from the mean state durations, and 
three from the frequencies of occurrence of the four states, which we estimate from information in KLIO. 
The forbidden transitions are represented by four zero elements, so only one additional equation is needed. 
We estimated mean durations (Ti,r2,T3,T4) ^ (1-01, 16.6, 1.01,8.9) (in period units) and frequencies of 
occurrence for the four states, / = (/i, /2, fs, fi) ~ (0.01, 0.655, 0.013, 0.321), from Figure 3 of KLIO. 

The full transition matrix was determined by searching over a grid of values for the off-diagonal elements. 



ql3, q23, q31, and q41, while holding the diagonal elements qi; 



1 , 4 fixed to values consistent with the 



mean durations and setting qi2 = q2i = 934 = ^43 = 0. For each trial transition matrix Q we determined 
the state probability vector P from Q (as an approximation to Q°°; see Sectional). By minimizing the 
mean-square difference |P — /p we obtained a (non- unique) solution. 



Qim 



(ni,n2,bi,&2) 



/ 0.01 0.886 0.104 \ 

0.940 0.006 0.054 

0.469 0.521 0.01 

\ 0.012 0.100 0.887 / 



(9) 



1989 1 in which the underlying physical states are manifested indirectly. We note that such models have been used to identify 



state changes in a pulsar ( Karastergiou et al. 12011 1 
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The values of 0.01 for the inetastabihties of short nulls and short bursts yield durations of order one time 
step, i.e. T = 1/(1 — 0.01) ~ 1. Higher time resolution data can sharpen these values. Figure [2] (left) 
shows the Markov distribution that superposes the two burst states (properly weighted according to their 
frequencies of occurrence and durations) overplotted on the histogram of burst durations from KLIO. The 
Markov distribution is a good representation of the salient features of the histogram. There may be an 
excess of observed long bursts compared to the observed counts, but they are small in number; some of 
these may result from short nulls that are missed in some of the data (i.e. false non-transitions). 

As implied in our merging of the two burst states, the four-state description for B1944+17 can be reduced 
to a three-state description if short and long nulls are combined into an S^ = {n, 61,62} niodel or bursts 
combined into an S' = {ni,n2,6} model. Similar 3-state descriptions can account for the dual nulling 



patterns seen from pulsars B0809-I-74 (van Leeuwen et al. 2002 Gajjar et al. 2012); B0818— 13 (Janssen 



& van Leeuwen 2004); J0828-3417, J0941-39, and J1107-5907 (Burke-Spolaor et al. 2012); B0835-41 



B2034-M9, B2021-h51, and B2319-h60 ( [Gajjar et alT||20l2l >. The object B0826-34 shows a null state with 



occasional single pulses and a burst state with rare single-pulse nulls (Esamdin et al. 2012) for which 



the average pulse shapes in these states differ significantly. This object is also amenable to a Markov 
description but a quantitative description of state durations and fractions is needed. 



4.3. Subpulse Drift Modes 

In some pulsars, there are well-organized motions of subpulses across the pulse-phase window defined by 
averages of a large number of pulses. Two or three subpulses typically appear in a single pulse period 
with separation P2, usually expressed in time units. In a sequence of pulse periods, bands of drifting 
subpulses have separations P3 at fixed pulse phase, usually expressed in units of the spin period P. The 
drift rate is then P2/P3 or, in physical units, (p = P2/P^Pi cycles s~^ when P2 and P are expressed in 
time units and P3 in spin periods. A common interpretation is that drifts correspond to E x B motions 
of multiple particle and radiation beams (a beam 'carousel') around a magnetic axis or some other related 



axis dRuderman k Sutherland|1975t[Deshpande fc Rankin|1999||Gil et al.|2006HLi et al.|2012b[ ). The total 
time needed for a beam to circulate around the axis is P4 = 1/4>P = PP3/P2 in units of spin period^ 
While P2 and P3 are visually recognizable in pulse sequences for some pulsars, the circulation time P4 



is often inferred from power spectra of pulse intensities at a single pulse phase (Backer 1973; Wright Sz 



Fowler|[l981b Deshpande Sz Rankin||1999 Weltevrede et al.||2006 ). For some objects discussed below, the 
drift pattern can have variable P3 and P4 but with near constancy of P2- 

When in the 'on' state, B1944-I-17 shows additional complexity in the form of four different emission 
modes (Deich et al. 1986| KLIO) that are distinguished by different subpulse drift rates or average pulse 
shape. Combined with nulling described above, as many as ten states may be needed for a full description, 
but additional observational classification is needed for any detailed modeling. 



Two pulsars that show three distinct drift rates along with nulls are B0031— 07 (Huguenin et al. 1970 



Wright fc Fowler||1981at [Vivekanand k Joshi||1997D and B2319+60 ( [Wright fc Fowler||1981b| ). The states 
are labelled A, B, and C in order of increasing drift rateQ For B0031— 07, the drift rates have ratios 
1:1.9:2.8 and observed sequences include B only and hybrid AB and BC bursts. There are no A only 



Pi is sometimes denoted P3 (e.g. Ruderman fc Sutherland|l975 Rankin fc Stappers[[2008l, which can be confused with 



P3 and is therefore not used here. 



*For B2319+60, we denote the 'abnormal' (ABN) mode defined by Wright & Fowler (1981b I as 'C 



13 



or C only bursts, or AC, BA, CA, or CB combinations. A four-state model S = {n, A, B, C} will have 
a transition matrix with q^c = QAn = QAC = Qba = QCA = QCB = and six remaining unique values 
need to be determined from measurements. Three can be obtained from the frequencies of occurrence 
/ = ifn, Ia, Ib, fc) ~ (0.45,0.086,0.45,0.014) and the four diagonal elements from the as-yet unknown 
mean durations of each state. Future measurements can allow determination of the full matrix with 
one redundant estimate that can be used as a check. Our discussion comprises an existence proof for a 
Markov model because, even without knowing the diagonal elements of Q, it provides the mechanism for 
'forbidden' transitions as well as the allowed transitions. 

Very similarly, B2319-I-60 has A,B, and C states that display different drift rates and pulse shapes along 
with nulls occurring with / = {fn, Ia, fs, fc) ~ (0.3,0.45,0.15,0.1). However, there are more allowed 
sequences than for B0031— 07, so eight unique elements of the transition matrix need to be determined 
after setting qac = Qba = QcA = QCB = 0. Seven can be estimated from the mean state durations 
(currently not in the literature) and /. The eighth can be obtained by estimating the number of AB or 
BC transitions. 

A remarkable feature of state transitions in both B0031— 07 and B2319-I-60 is that subpulse drift states 
occur only in order of increasing drift rate. In the carousel model, the drift rate is related to the local 
plasma drift velocity cE x B/B^ that vanishes when the local charge density equals the Goldreich-Julian 
density. An increasing drift rate therefore suggests depletion of the charge density over the course of a 
burst that eventually terminates in a null. In a nonlinear model, the systems in these two pulsars seem 
to be running 'downhill' in an effective potential that is somehow reset during and perhaps because of a 
null. The A, B, and C states would have energy minima such that it is much more likely to have AB, 
BC, and AC transitions than BA, CB, or CA transitions. 



4.4. Mode Changes with Long Time Constants 

Two pulsars show pulse-shape modes that have finite time constants for the evolution of the drift rate in 
one of the modes. The exemplar drifting-subpulse object B0809+74 has 1.4% null pulses that interrupt 



bursts (Lyne & Ashworth 1983; van Leeuwen et al. 2002). Like other objects, the durations of nulls and 



bursts have Markov distributions, but during bursts the drift rate decreases at the onset of a null and 
relaxes exponentially to a larger asymptotic value over a time that scales with the null duration but is 



typically tens of seconds (Lyne Sz Ashworth 1983). 



The B and Q burst modes from B0943-I-10 ( Rankin &: Stappers|2008 ) can be accommodated by a two-state 
Markov process but with some complications. This pulsar displays correlated X-ray and radio states that. 



along with torque variations seen in other pulsars, link radio and high-energy phenomena (Hermsen et 



aL]|2013 ). The two burst-only radio states show a quasi-period P4 that is constant in the weaker Q mode 
but evolves exponentially with a 73 min time constant in the B mode. As P4 increases the pulse shape 
changes progressively on the same time scale. The exponential increase in P4 (equivalent to a decrease 
in drift rate) suggests in the carousel model that, unlike for B0031— 07 and B2319-I-60 discussed above, 
the E X B drift velocity decreases in B-mode sequences, implying that the polar-cap acceleration region 



trends toward a force-free state. In the ionic evolution model of Jones (2011), the vertical potential drop 



in the polar-cap accelerator evolves as high-Z ions are increasingly ionized and disassociated. It is not 
clear in either of these interpretations why the state changes are not more periodic rather than showing 
a wide range of durations, but a nonlinear system driven stochastically can account for these. 
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5. Long-term Intermittency and Stochastic Resonance 

Three objects (J1832+0029, J1841-0500, and B1931+24) show distinct on and off emission states that 
are accompanied by switching between two values of the spindown rate v. Thus far there are no detailed 
descriptions of pulse variations in the on state because intensity and timing results have been based on 
average profiles rather than single pulses. 

The best studied object, B1931+24, shows no short-term nulls (minutes or less) when the pulsar is in 



the long-term 'on' state (Kramer et al. 2006). The durations of ons and offs have mean values of 8 ± 4 d 



and 22 lb 7 d and large maximum to minimum ratios, 19:1 and 10:1, respectively (Young et al. 2013). 



These ratios are similar to those of short-term nulls and bursts for other objects and, along with the lack 
of correlation between the durations of contiguous on-state and off-states, are consistent with a Markov 
description. 
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Fig. 6. — Left: Histograms of on-and-off state durations for B1931+24, after Figure 2 of Young et al. 
(2013). Vertical dashed lines indicate one-half and a full period of Pj = 38 it 4 d identified in a 13 y 



data set (Young et al. 2013). Right: Time series of state changes and histograms or state durations for 
a stochastic resonance model of the form of Equation [3] with q\\ = 0.985,^22 = 0.83, Ai = 2.6, A2 = 1.6, 
and Pf = 38 d. 



However, the histograms of state durations are inconsistent with the monotonic form expected from a 



Markov process. As shown in Figure ro^ (left panel), based on Figure 2 of Young et al. (2013) (see also 



Figure Ic of Kramer et al. (2006)), the mean state durations are about equal to their modes, and there 
is a paucity of short durations. The overall small number of on-off cycles in the 13-yr data set combined 



with sampling incompleteness (from allowing observation gaps as large as 5 d in the analysis. Young et 



al. (2013)) may contribute to the on-state histogram shape. Conversely, the deficit of off-state durations 
up to ~ 15 d is not plausibly accounted for this way. 

Another difference from fast intermittency and Markov behavior is an apparent underlying period for 
the intermittency. Pi = 38 it 5 d, in the 13 yr data span (Young et al. 2013) that is longer than any 



of the individual on-and-off durations and is about 4.8 times the mean on-duration and 1.7 times the 
mean off-duration. The quasi-periodicity may be interpreted in two ways. First, there may be a strictly 
periodic — but hidden — process that is blurred by mediating processes that lead to observable quantities. 
Alternatively, the 38 d period may be the mean of a process that is intrinsically quasi-periodic. We explore 
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the first of tliese options because it is more specific. 

A stocliastic-resonance model is an acceptable statistical description of observations of B1931+24 if we 
adopt a transition matrix (Equation^ with a sinusoidal forcing function d{t) = sm{27rt/ Pf + </>). In this 
interpretation the different mean on-and-off durations Ton,T'off signify an asymmetric effective potential 
and their sum, Ton + Toff = 31 zb 8 d, is nominally smaller than the observed mean period P, = 38 it 5 d 



reported by Young et al. (2013). Schmitt triggers show similar behavior when the underlying potential 



has asymmetric minima ( Marchesoni, Apostolico &: Santucci||1999 ). 

The right-hand panel of Figure [6] shows state-duration histograms for a stochastic-resonance model that is 
consistent with the observed histograms in the left-hand panel. The simulations include a strictly periodic 
forcing function with Pj = Pi. The centroids of the dominant features in the simulated histograms straddle 
the half period Pf/2 = 19 d, as is common in SR models. There is an additional feature in the simulated 
off-state histogram for durations < 7 d that may be absent in the actual data because of the incomplete 



sampling of short-duration states (Young et al. 2013) and because the total number of counts is smaller 



than in the simulated histogram. A more precise test of stochastic resonance in B1931+24 would require 
longer, continuous data sets than are currently available. Data sets that include single-pulses would be 
valuable for characterizing any very short-duration on or off states. 

The other two objects that show long-term intermittency and discrete values of z>, J1832+0029 and J1841- 
0500, have state durations that are longer than those in B1931+24 and have not been sampled adequately 
to test whether a stochastic resonance model might apply to them as well. 



6. Timing Variations from a Tvi^o-State Markov Process 

Objects that are intermittent on long time scales like B1931+24 show a distinct value of spindown rate 
in each state. Switching between discrete values of z> can be viewed as pulse-like changes in z> (i.e. equal 
upward and downward transitions) that combine with a constant value of z>, where we assume that the 
second derivative i> from the spindown torque is negligible. Viewed collectively as a two-state Markov 
process, switching between spindown rates produces a random walk in spin frequency and an integrated 
random walk for pulse phase. Stochastic timing variations occur because individual states last for times 
that vary widely about their means, as discussed earlier. We integrate the fluctuating part of z>(t) twice to 
get the spin phase perturbation and calculate its variance over a data span of length T^^ T . Using the 
'reduced' duration, T = TiT2/{Ti + T2), the rms timing variation expressed in time units can be written 
as 



cTtMTd) = C^^MTd) « 74 fisP\u^i5\ tJ/'tJ/'^ 



50 
ii 



Vfi{i-fi) 

1/2 



(10) 



In this equation, we have used T = 10 Ts s, Td is in years, and z> = z>_i5l0^ s~ . The factor C2 ~ 15.5 
corrects for the fluctuations removed by subtraction of the second-order polynomial fit for the mean spin 
rate and z>. All else being equal, the rms timing error is a maximum for equal fractions of low and high 
states, fi = fh = ^ — fi = 1/2. The scaling, at{Td) on TJ , is generic to random walks in spin frequency 
( |Groth|[T975l [Cordes fc Dowslp85| [Shannon fc Cordes|[2010l ) . 

For B1931+24, individual transitions can be identified and removed from the data to reduce the rms 
timing variation ( Kramer et al.|2006 ). However for other objects, individual states may not be identifiable 
because they are too small or are mixed with other spin variations and consequently will appear as spin 
noise with nonstationary variance. 
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7. Medium-term Intermittency and Timing Variations in B0823+26 

B0823+26 (P = 0.53 s and z> = 6 x 10~^^ s~^) displays minute-duration nulls along with those of hours 
and longer. Any emission during nulls is less than 1% of the mean on-state flux density. The distribution 
of on-state durations shows both a narrow and a broad component with times scales of minutes and days, 



respectively (Young et al. 2012). The durations of particular nulls and bursts may be mis-estimated in 



some cases due to data windowing effects that cause state changes to be missed, as noted by Young et al 



(2012), but the presence of both long and short nulls is secure. 



Torque variations are evidently more complex in this object than in long-term intermittent objects like 



B1931-I-24. In a 153-d data set. Young et al. (2012) saw neither cubic structure in the timing residuals 



above an rms residual of about 200 fis nor switching between discrete values of u like that seen in 
B1931+24. The upper bound on the amplitudes of such pulses is Au/u < 0.06 ( Young et al.||2012 ) based 
on detailed modeling of the time series. What are seen in a longer data set, however, are step-function-like 
changes in z> (rather than pulses) with amplitudes as large as |Az>/z>| = 0.02 but occurring at a rate of 
only ~ 1.2 yr~^ ( Cordes Sz Downs|l985 ), much less frequent than the state changes discussed here. 

We can place an additional statistical limit on the amplitudes of pulses in z>. Equation [To| with pi = 0.2, 
Ti = 0.26 d and T2 = 1.4 d ( Young et al.||2012 > predicts (Tt,2(153 d) « 610 /xs \du/u\. Comparison with the 
upper bound of 80 fis on the rms residual (c.f. Figure 8 and discussion in Young et al. (2012)) implies 



|5z>/z>| < 0.13, a less stringent constraint than \dO/i^\ ^ 0.06 obtained by (Young et al. 2012) from direct 



fitting. Using the smaller upper bound, we estimate that for a 13.6-year data set like that analyzed by 
Cordes Sz Downs (1985), state changes would produce no more than cr(^2(13.6 yr) « 6.7 ms, much smaller 

100 ms over this data span. Spin noise in B0823+26 is evidently dominated 

much steeper than that expected from pulses in 



than the measured value of ■ 

5/2 
by step functions in z> that yield a scaling (Tt(T^) oaT, , 



z>. The origins of the step functions in z> for B0823+26 appear distinct from the process that causes large 
pulses in z> for B1931+24 ( [Kramer et al]|2006 ). 

Prom the timing analysis on B0823+26 we conclude the following. First, switching between nulls and 
bursts on time scales of days evidently is not associated with large jumps between binary values of 0, the 



same conclusion as reached by Young et al. (2012). Second, the large amount of spin noise appears to 



have other causes than those that are responsible for the intensity jumps. The difficulty in identifying 
discrete values for z> in B0823+26 may result from a two-state model being too simplistic. Any correlated 
torque events might also have different shapes than pulse-like forms. If the torque does not return to its 
original value (as it does in B1931+24 and other intermittent pulsars), the spin noise will appear more 
like that observed. Multiple states are also implied by the presence of nulls and bursts that each have 
short and long-duration components. A timing analysis of multiple state models is beyond the scope of 
this paper but is worthy of future analysis. 



8. The Physics of Discrete States and Triggering Between States 

The aggregated observational results presented in this paper demonstrate that diverse phenomena result 
from the existence of metastable states in the dynamics of pulsar magnetospheres. A physical understand- 
ing of the phenomena leads to two primary questions: (1) What is the underlying magnetospheric physics 
of the observed states? and (2) What causes transitions between states? We discuss these questions in 
the context of standard ideas about pulsar magnetospheres. 

Some characteristic light-travel times associated with neutron stars (NSs) and their magnetospheres are 
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useful for the following discussion. The NS radius is /?=„ ~ 30 /zs and the radius of the PC is rpc = 
-R*(i?*/ric)^'^ ~ 0.5P^^'^ fis for spin periods P in seconds, where the light cylinder radius ri^ = cP/2tt ~ 
0.16P s. In models with an acceleration region just above the PC, the height of the region h < rpc 
and for objects with sustained pair cascades, h ~ 0.3 /is. Radio intensity variations are seen on these 
and smaller time scales. Indeed, a priori one might expect all observed time scales to be no larger. 
However, the variations discussed in this paper concern macroscopic scales that exceed the characteristic 
times by many orders of magnitude (seconds to years) but are still much smaller than the spindown 
time Tg = P/P ~ W^'^PPZi^ yr on which the magnetosphere evolves (where P = 10~^^P_i5 s s~^ is 
the period derivative). Thermal time scales enter the picture when thermionic emission plays a role and 
the temperature of a hot PC is time variable. Also, in some situations, pair cascades are temperature 
dependent ( Hibschman fc: Arons|[200T Medin & Lai 2007). The heating and cooling times of the PC are 
estimated to be < 1 s and thermostatic regulation from pair production may result in variations on some 



multiple of this time scale (Cheng Sz Ruderman 1980). Surface physics introduces the proton drift time 
through the atmosphere Tp ~ 0.1-1 s and the 'excavation' time for removing one radiation length's amount 
of material from the PC, th ^ 2.1 x 10^ s {P/ZB12) ( |Jones||201l||2012[ ). Orbital periods around a 1.4 Mq 
NS are Porb ~ 127 sP'^'^(r/ric)^'^ using r\c as a fiducial radius. An object with a 38 d orbital period at 
r ~ 0.3 AU is well outside the gravitational tidal disruption radius, rtg = (3MNs/27r/9)-'^'^ ~ lO^^p"-*^'^ cm. 
Small objects inside the tidal radius can exist owing to tensile forces. 

Radio emission requires relativistically beamed, coherent radiation from a particle flow that is collimated 
in a narrow, magnetic field-line bundle. Coherence is produced by an instability that requires counter- 
streaming particle flows and has a maximum growth rate at an altitude-dependent plasma frequency. Also 
required is a favorable viewing geometry that depends on the angular width of the radiation beam, which 
varies with altitude, and on the angle between the beam and the observer's direction. The observed 
discrete states therefore require changes in one or more of these elements to alter the intensity, pulse 
shape, and spectrum of the radiation. That discrete states are seen in the dominant energy loss channels 
— viz. spindown rates v in the long-term intermittent pulsars like B1931+24 and X-ray emission from 
B0943+10 — implies that significant, discrete changes occur in the magnetosphere at large. 



8.1. Metastable States 

Several gross features of pulsar magnetospheres could provide binary states, including on-and-off mod- 
ulation of e production, cycling of the relative contributions of ions and pairs to the current flow, or 
modulations of the total current itself between two extremes, such as a vacuum state (no current) and a 
force-free state (e.g. Li et al. 2012b|c ). Another dyadic pair could be a dead, 'electrosphere' state (e.g. 
Michel|1980t [Krause-Polstorff k Michel|1985[ [Smith et al.|2001t [Petri et al.|2002p vs. a standard, dynamic 
magnetosphere. Extreme switching between an electrosphere and a magnetosphere might be associated 
with the more extreme RRATs and long-term intermittent objects (e.g. 



Michel 2010). 



Observed state changes so far are exclusively the realm of older, canonical pulsars with periods > 0.3 s 
and magnetic fields B ~ 10^^ G. Indeed, there appears to be a greater propensity for state changes in 
pulsars that are longward of ~ 0.3 s in the P-P diagram (Figure 1, Cordes &: Shannon 2008). However, 
it is conceivable that shorter period objects display state changes that have been missed because of 
observational selection, particularly MSPs for which there are few observations of single pulses. Whether 
any state changes are expected from MSPs depends on the nature of the states themselves and on how 
state changes are triggered. If one of the state dyads is associated with a threshold for pair production. 
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state transitions will be expected if the pair-production death line is near the MSP population in the P-P 
diagram. Published death lines do not provide clarity on this because they depend on the types of current 



flow (SCL or not), on the photon-emission process (curvature vs. inverse-Compton) (Hibschman & Arons 
2001t [Harding et al]|2002t [Wang k Hirotani||2011t [Harding fc Muslimov||2011t [Medin fc LaT[|2007[ ) , and on 



any offset of PCs from the dipole axis caused by current-driven sweepback of the magnetic field (Harding 



&: Muslimov||20lT ) or by an offset dipole (e.g. Hibschman &: Arons[[200T ) . If alternative photon-processes 
can drive pair cascades for different regions in the P-P diagram, the propensity for state changes may be 
associated with just one of the processes, such as non-resonant inverse Compton radiation vs. curvature 
radiation. 



8.2. State Triggering 

A model is proposed in which the dynamics of current flows are responsible for both defining discrete states 
and for causing transitions between them. In this picture, spindown rates and X-ray states, along with 
radio nulling, drifts, and mode changes, are all collateral effects. Discrete states are tied to the physics 
of the magnetic polar cap (PC) while switching involves both local and global effects. The suitability 
of simple Markov processes for describing the distributions of state durations in most cases suggests 
that state switching is largely stochastic. This is most easily understood as self-driven stochasticity of 
the acceleration region. The natural time scale for such variations < 1 ^s. However, the stochastic 
resonance model for the quasi-periodicity in B1931-I-24 suggests the action of a forcing function with a 
long period ~ 38 days that is easiest to understand in terms of modulation of the return current from an 
equatorial disk outside the light cylinder. Equatorial disks are a common feature in recent discussions of 
magnetospheres and in numerical simulations of pulsar magnetospheres ( Krause-Polstorff &: Michel|[1985 



Contopoulos 2005 Timokhin 2010 Kalapotharakos et al. 2012 Li et al. 2012b|c) and have long been 



argued for by Michel (1980). Dramatic changes of magnetospheric structure by pulsar-disk interactions 



were a conclusion of Barsukov et al. (2009). Recent magnetosphere models incorporate a 'Y' point where 



current from a disk flows along the separatrix to both magnetic PCs (e.g. Contopoulos 2005 Timokhin 



2006). This feature may provide the most direct feedback from a pulsar's environment to its PCs. 



8.3. Circuit Analogs 

The NS surface and atmosphere at the PC and the acceleration region above have been likened to the 



idealized unidirectional current flow of a diode (e.g. Arons & Scharlemann 1979 Thompson 2008). The 
analogy has been extended to include resistivity, capacitance, and inductance in linear circuits (e.g. 



Shibata 1991; Jessner et al. 2001 Xu et al. 2006) that, in principle, could show oscillations with well- 



defined time scales. However, nonlinear circuits can have metastable states with broad, Markov-type 
durations when state switching occurs from noise-like voltage variations. A nonlinear, real-world diode 
with a forward bias voltage, a back-biased capacitance, and a finite recovery time to a reversal from forward 



to backward bias can show universal chaotic behavior when driven sinusoidally (Rollins &; Hunt 1982). 



Nonlinear circuits are amenable to a Markov description when switching between states is rapid compared 
to the residence time in any state and when 'microscopic' dynamics within a state can be ignored compared 
to the step-size of the Markov chain (e.g. Becker &: Rein ten Wolde|2012 ). An asymmetric Schmitt trigger, 
which features hysteresis and nonlinearity, has unequal switching thresholds for positive and negative 
inputs. It shows Markov variations for a random-noise input but can show stochastic resonance when 
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driven sinusoidally, behavior that is generic to any asymmetric bistable system (Marchesoni, Apostohco 



k Santucci 1999). 



NS and their magnetospheres differ substantially from laboratory circuits because voltage drops and 
current flows along field lines are not controlled but instead are dynamical quantities, providing a richer 
variety of possible phenomena. 

Pulsar magnetospheres strive to be force free {E ■ B = 0) hy filling themselves with space charge at 
the Goldreich-Julian (GJ) density pcj ~ —^ ■ B/lnc, where ii is the spin vector. Standard models for 
pulsars include a corotating magnetosphere with p = pQj and an open field line (OFL) region where 
charge is lost through a current density J ^ cpcj, but p = pcj almost everywhere except in acceleration 
regions where E ■ B ^ 0. The potential drop sustains the current flow that contributes to the spindown 
torque and, along with secondary particles, produces radio to gamma-ray emission. The OFL region and 
magnetic PC are defined self-consistently from the combined ambient stellar magnetic field and current- 
driven field; the latter, in turn, depends on boundary conditions imposed by the global structure of the 



magnetosphere (e.g. Kalapotharakos et al. 2012). 

Recent simulations ( Kalapotharakos et al.|[2012 Li et al. 2012b|c ) demonstrate exquisitely how the spin- 
down rate depends on the magnetic flux threading the light cylinder and takes on limiting values for 
vacuum (p = 0) and force- free {p = pcj) conditions in the OFL region. Spindown rates in these limits 
differ by a factor of four for inclination angles of 45°, easily bracketing the 2.5:1 range seen in long-term 
intermittent pulsars. For active, steady pulsars the magnetosphere appears to find a configuration that 
lies between the force-free and vacuum solutions. However, there is no obvious mechanism for defining 
dual (let alone multiple) states in what appears to be a continuum for the resistivity, which is a free 
parameter in determining the type of magnetosphere ( Li et aiy]|2012b|c ). 

We attribute discrete current states to changes in the relative contributions from ions, protons, electrons, 
and e pairs. This situation applies to an fi • S < geometry for the OFL region at the NS surface, 



the situation considered by Ruderman Sz Sutherland (1975) and Jones (2011 2012 2013), that requires 



positive charges above the PC. The opposite case with Q ■ B > has comparatively uncomplicated 
electron flows that may not allow pair production and radio emission. However, if it does, state changes 
are not expected. The ft ■ B < and Q ■ B > cases may then account for the stark differences between 
pulsars having similar P and P. Positive charges include ions with a wide range of charge to mass ratios 



resulting from particle showers onto the polar cap (Jones 2012) along with protons and any positrons 



produced in pair avalanches. One form of current self-regulation involves thermionic emission from a PC 
that is kept hot by backflowing electrons from e^ pairs (Cheng & Ruderman 1980). Small changes in 



temperature can change the ion current density Jj exponentially if the current flow is not space-charge 
limited (SCL). In this case, cessation of e production allows the PC to cool and Ji will decrease. Unless 
another charge source can compensate, the voltage drop will increase while the spindown rate decreases. 
Coherent radiation that requires a substantial pair plasma will also shut off. With SCL flow, however, 
cessation of pair production will cause a radio null but without a large decrease in torque on the star. 



8.4. Return Currents and Stellar Charge 

Current outflow must be matched by a return current, at least on average, if the NS is not to become 
charged sufficiently to halt current outflow from the polar cap. While a disaster for explaining steady 
pulsar radiation, episodic charging and discharging is a plausible mechanism for intermittency because it 
could take little time (about one spin period) to switch between charge states. In some models the return 
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current is from separated e^ pairs created within the magnetosphere, while in others it originates from an 
equatorial disk exterior to the magnetosphere. If a mismatch develops between outflow and inflow, it will 
self regulate if pair production is robust as in high-field, fastly spinning pulsars. However, longer-period 
or weak-field pulsars that operate near the threshold for pair production, will have a higher propensity 
for state changes. It is therefore notable that fast states (classical nulling, drifting and mode changes), 
long-term intermittency, discrete torques, and X-ray states have been seen only in pulsars with P > 0.3 s 
and are common for P > 1 s. Furthermore, nulling is most visible in objects with small angles between the 
spin axis and magnetic moment ( Biggs||l992 Cordes fc Shannoii||2008 ) . Given both of these conditions, it 
is plausible that currents originating from an equatorial disk can drive the system into and out of a state 
of dormancy, at least as far as radio emission is concerned. 

Return currents from an exterior disk, if highly episodic, will modulate activity at the PC as well as alter 
the structure of the magnetosphere ( Barsukov et al.|2009 ). Depending on their constituency, they can also 
'dope' the NS surface with ions that differ from those that drift to the top of the thin atmosphere from 



below (Jones 2011), which are fractionated in charge-to-mass ratio by backflow electrons (Jones 2012). 



In laboratory diodes with non-SCL current flow, the current density depends on the type of material 



as well as on the voltage drop (Goncalves, Barroso, & Sandonato 2004). A similar effect may occur in 



pulsars where the charge carriers could be a state-dependent mixture of surface and backflow ions. Like 
any circuit, the current flow in the magnetosphere can halt altogether. With sufficient net charge, the 
combined electrostatic and induced potential drops can shut off the conventional PC accelerator by reverse 



biasing the system (Petri et al. 2002). 



These considerations suggest that if return currents originate from or are affected by a disk external to 
the light cylinder, the electrodynamics of the disk itself will be involved in state-change phenomena. Disk 



instabilities combined with perturbations from contaminating material (e.g. asteroids; Cordes & Shannon 



2008) may introduce a wide range of time scales, including orbital time scales. Some of these may be 



stochastic and others periodic. In addition, magnetic reconnection within the disk from magnetic fields 
with oppositive polarities above and below the disk can inject relativistic particles into the magnetosphere 



(e.g. Contopoulos 2005; Komissarov 2006 Timokhin 2006). 



8.5. Discrete Drift States 



We associate fast state changes with the acceleration region just above the NS surface and atmosphere e.g. 



(Ruderman &: Sutherland 1975). Fast nulling and pulse-shape modes suggest alteration of the particle 



composition in current density and the number of secondary e pairs. Preferred drift rates may be 
associated with changing numbers of sub-beams in a rotating carousel. The number of sub-beams is 
likely involved in determining the local charge difference Ap = p — pcj and thus would affect the E x B 
drift velocity. The number of carousel beams ~ 20 in B0943-I-10 (Deshpande & Rankin 1999) while 



another object, B1822-09, appears to involve only two subpulse beams in B mode and three in Q mode 



(Latham et al. 2012). Sub-beam creation and annihilation would be an attractive process for producing 



discrete values of subpulse drifts and other state attributes. However, two pulsars mentioned earlier 
(B0031— 07 and B2319-I-60) have approximately constant P2 in the three observed drift modes in which 
there is an increase in P3 and a decrease in the drift rate P2/ Pz- The spatial separation of sub-beams is 
directly related to P2 suggesting, at least for these two objects, a constant number of sub-beams because 
Nsb oc Pa/ P2- In other pulsars the variability of P2 is not well constrained. An alternative to the rotating 



carousel model is the diocotron instability that has been suggested as a cause of subpulse drifts (Fung et 
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ar]|2006 ), which results from differential angular rotation across the OFL region. Magnetospheric models 
that incorporate this instability are not well-enough developed to provide a basis for further discussion 
in this paper. 



8.6. Long-lived States, Energetics, and MSPs 

The range of time scales seen in pulsar intermittency is large so it is natural to ask what other time 
scales are relevant and what the energetics are. Given that observed time scales have little relationship 
to characteristic times of NS and their magnetospheres and that processes exterior to the magnetosphere 
may be involved, objects with much longer state durations (e.g. years, decades) plausibly exist that will 
be discovered in comprehensive radio surveys that have a large product of time coverage and field-of-view. 
(e.g. [Macquart et al.||2010[ [Murphy et aI]|2013D . 

As for energetics, it is well known that pulsar radio luminosities are only small fractions (~ 10~^ to 10~^) 
of the spin-down loss rates E but can be comparable to particle energy fluxes. It therefore follows that 
sporadic on-state emission can be seen over distances that are limited by the amount of energy available to 
a burst of particle acceleration. If magnetospheres make occasional transitions from a near vacuum state 
to a force-free state and back, the available energy can be much greater than for typical radio pulsars and 
the resulting pulses may be visible over extragalactic distances. Such pulses might be related to any that 
may be seen with dispersion measures too large to be accounted for by sources embedded in the Galactic 



disk (e.g. Keane et al. 2011). 



9. Summary and Conclusions 

In this paper it has been shown that metastable states in pulsar phenomena are common and that while 
many cases involve binary states, some pulsars display multiple states. Pulsars are well known to show 
pulse profiles that are peculiar to each object but are consistent over tens of years or more. Similarly, 
the presence and properties of state switching also appear to have stationary statistics. These properties 
include the mean durations and the long-term frequencies of occurrence of particular states, including the 
nulling fraction often used to characterize objects that show intensity nulls. Collectively, the phenemona 
suggest that states are defined by long-lasting features of neutron stars and their magnetospheres. In the 
Markov interpretation presented in this paper, the transition matrix appears to be epoch independent in 
most objects other than B1931+24. 

We have also shown that quasi-periodic switching on time scales of weeks for B1931+24 is consistent with 
a Markov system that displays stochastic resonance. Stochastic resonance is seen in systems that are 
driven by both stochastic fluctuations and a forcing function. Most studies assume the forcing function is 
sinusoidal, but arbitrary deterministic or stochastic forcings are also generally possible. For B1931+24, 
the observed quasi-periodicity allows either a deterministic or quasi-periodic forcing function. It is not 
obvious how to discriminate between the two observationally. One possibility is to monitor this object 
(and other long-term intermittent pulsars) with high-cadence observations so that switching statistics can 
be better compared with specific models. 

To date, only a small fraction of pulsars have been studied with sufficient S/N to allow quantitative eval- 
uation of a Markov (or other) model. Even for many of the well-studied objects, pulse-to-pulse variations 
cause overlap of on-and-off-state intensities, leading to false positives from algorithms that identify state 
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changes. Future work with now-available wideband receiver systems on existing telescopes and eventually 
with new array telescopes (ASKAP, MeerKAT, and the SKA) can improve the discrimination between 
states as well as expand the sample of objects that can be studied in this way. 

Of equal interest are objects that do not show (or have not yet shown) any kind of state change. Pulsars 
with periods < 0.3 s, including MSPs, are among these. Such pulsars may have such robust pair produc- 
tion that they are not susceptible to state changes because return currents are self-generated inside the 
magnetosphere rather than originating from an external disk. In addition, there are objects that do not 
show state changes yet are in the same region of the P-P diagram as pulsars that do. It is important 
to establish more firmly the occurrence or non-occurence of metastable states in these pulsars. If pulsars 
that are otherwise similar in P and P show markedly different variability, additional factors need to be 
considered that may be peculiar to individual objects. Candidate factors include the angular offset of the 
magnetic axis from the spin axis, higher-order multiples of the magnetic field, surface composition, and 
environmental effects in or near the equatorial disk. Another possibility is that the charge polarity at the 
PC, given by the sign of fi • B, determines whether a pulsar displays metastable states. 

Future work that can illuminate some of the issues raised in this work include: 

1. Simultaneous observations over wide frequency ranges to test whether some null states represent 
redistribution of radio flux to a different frequency band; 

2. Long-term monitoring campaigns with single-pulse sampling to identify state changes that occur on 
both short and long time scales; 

3. Characterization of state changes as Markov processes in more pulsars and testing for long-term 
stationarity of the transition matrix; 

4. Detailed investigations of single pulses from MSPs to assess whether state changes occur that are 
qualitatively similar to those of canonical pulsars with '^ 10^^ G fields; 

5. Higher sensitivity observations that discriminate between null and burst states with minimal or no 
false transitions; these can also clarify whether there are true null states or simply low-radio-emission 
states; 

6. High cadence timing measurements that allow differentiation between state changes and other con- 
tributions to spin noise; 

7. Very low frequency observations to test whether there is coherent ion emission; and 

8. Additional simultaneous high-energy and radio observations to elucidate the nature of emission 
states and the locations of the emission regions for the different wavebands. 
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